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ABSTRACT 

Ionized gas discs in the nuclei of ellipticals have proven to be excellent tools for the 
£j>^ | determination of the central mass density in these galaxies. The recent discovery with 

the Hubble Space Telescope of small stellar discs embedded in the nuclei of a number 
of ellipticals and SOs might be of similar importance. We construct two-integral ax- 
^ \ isymmetric models for such systems. The models consist of a spheroidal bulge with a 

central density cusp, and a disc described by a strongly flattened exponential spheroid. 
I/*-) ' We use the Hunter & Qian (1993) method to calculate the even part of the phase-space 

distribution function (DF), and specify the odd part by means of a simple parameter- 
ization. We consider both local stability against axisymmetric perturbations, as well 
as global stability against bar forming modes, and find that our models are stable as 
, long as the discs are not too flat and/or compact. The margin of stability is derived 

as a function of disc scalelength and central surface density. Its location agrees well 
with the observed values of these disc parameters. This suggests that discs build up 
their mass until they become marginally stable. 

We investigate the photometric as well as the kinematic signatures of nuclear discs, 
, including their velocity profiles (VPs), and study the influence of seeing convolution. 

In particular, we study to what extent these kinematic signatures can be used to 
determine the central density of the galaxy, and to test for the presence of massive 
black holes. We consider nuclear discs that are either dynamically coupled to or 
Q ■ decoupled from the host elliptical, including counter-rotating discs. The latter are 

Jjj \ models for the counter-rotating cores observed in a number of galaxies, which are 

c/3 . often found to exhibit discy isophotes in the central region. The counter-rotation is 

only detectable when the disc light contributes significantly to the central velocity 
profiles. We find that in this case the observed velocity dispersion will show a central 
decrease. 

' The rotation curve of a nuclear disc gives an excellent measure of the central mass- 

to-light ratio whenever the VPs clearly reveal the narrow, rapidly rotating component 
associated with the nuclear disc. Steep cusps and seeing convolution both result in 
central VPs that are dominated by the bulge light, and these VPs barely show the 
presence of the nuclear disc, impeding measurements of the central rotation velocities 
of the disc stars. However, if a massive BH is present, the disc component of the VP 
can be seen in the wing of the bulge part, and measurements of its mean rotation 
provide a clear signature of the presence of the BH. This signature is insensitive to 
the uncertainties in the velocity anisotropy, which often lead to ambiguity in the 
interpretation of a central rise in velocity dispersion as due to a central BH. 

Key words: stellar dynamics - galaxies: kinematics and dynamics - galaxies: ellip- 
ticals - galaxies: structure - galaxies: nuclei - line: profiles 
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1 INTRODUCTION 

Ever since the realization that elliptical galaxies are not 
smooth structureless systems of old stars supported by rota- 
tion, numerous studies have revealed the increasingly com- 
plex nature of these objects (de Zeeuw & Franx 1991). In 
addition to the bright, pressure supported triaxial systems 
with boxy isophotes, there are the discy ellipticals that are 
mainly supported by rotation. The central regions often 



exhibit complex structures. Counter-rotating cores, stel- 
lar cusps, and nuclear activity are all phenomena connected 
with the nuclei of elliptical galaxies (Kormendy & Richstone 
1995). In addition, the nuclear regions might conceal mas- 
sive black holes. 

Numerous studies over the last decade have concen- 
trated on the nature of the discy ellipticals (e.g., Capaccioli, 
Held, & Nieto 1987; Carter 1987; van den Bergh 1989; Ni- 
eto et al. 1991; Rix & White 1990, 1992; Scorza & Bender 
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1990). Scorza & Bender (1995) have shown that discy el- 
lipticals follow the same trend as observed for spirals and 
SOs that discs with smaller scale lengths have higher central 
surface brightness (e.g., Kent 1985). This strongly suggests 
that discy ellipticals are two-component systems that form 
a continuing sequence in the Hubble diagram from SOs to 
galaxies with smaller disc-to-bulge ratios (see also Michard 
1984; Capaccioli, Caon & Rampazzo 1990). 

The Hubble Space Telescope (HST) has revealed very 
small (scale length ~ 20pc), nuclear discs in a number of 
ellipticals (van den Bosch et al. 1994; Forbes 1994; Lauer 
ct al. 1995). Several of these galaxies also harbour a much 
larger kpc-scale disc which may have an inner cut-off ra- 
dius at 200 — 400pc. Such multiple component discs have 
been observed in other galaxies as well (e.g., Seifert 1990; 
Scorza & Bender 1995). For example, Burkhead (1986) and 
Kormendy (1988) argue that the huge disc of the Sombrero 
galaxy (NGC 4594) consists of two separate components. 
Emsellem et al. (1994) even find an indication for a third 
disc component in the Sombrero, namely a red nuclear disc 
that remains unresolved from the ground. 

In this paper we investigate the dynamical proper- 
ties and observables of nuclear discs. We construct self- 
consistent axisymmetric dynamical models of a spheroidal 
bulge with an embedded nuclear disc. In particular we in- 
vestigate whether the observable dynamical properties of 
nuclear discs can be used to discriminate whether or not 
these galaxies harbour a nuclear black hole (BH). Proving 
the presence of a BH in the nucleus of an elliptical galaxy 
is complicated by the fact that the three-dimensional ve- 
locity distribution of the stars in the central region can be 
strongly anisotropic. For example, the observed central in- 
crease of the velocity dispersion in M 87 can equally well 
be fitted by an isotropic model with a central BH as by 
a radially anisotropic model without any BH (Young et 
al. 1978; Sargent et al. 1978; Binney & Mamon 1982; van 
der Marel 1994). Discs, on the other hand, are cold compo- 
nents, strongly dominated by rotation. Furthermore, they 
are generally thin, so that the interpretation of the observed 
line-of-sight velocity distribution is less complicated than 
for the three-dimensional distribution of stars in the ellip- 
soidal component. Discs embedded in the nuclei of elliptical 
galaxies offer a promising means of detecting nuclear BHs. 
The recent discovery of a disc of ionized gas in the nucleus 
of M87 (Ford et al. 1994) has indeed greatly strengthened 
the case for a nuclear BH in this giant elliptical (Harms et 
al. 1994). Even more spectacular evidence for the presence 
of a nuclear BH is provided by the VLBA observations of 
emission of a disc of water masers with a Keplerian rotation 
curve at only 0.13 pc from the centre of NGC 4258 (Miyoshi 
et al. 1995). However, the interpretation of the dynamics 
of gas discs can be complicated, especially in AGNs where 
outflow, inflow, and turbulent phenomena are very likely to 
play a major role (see e.g., van den Bosch & van der Marel 
1995; Jaffe et al. 1996). Since the motion of stars is purely 
gravitational, stellar discs do not suffer from these effects, 
and their kinematics might offer valuable tests of the pres- 
ence of nuclear BHs. Recent absorption line spectroscopy 
with the Faint Object Spectrograph aboard the HST of the 
nucleus of NGC 3115, one of the galaxies with a bright nu- 
clear, stellar disc, has indeed provided an excellent case for 
the presence of a 2 x 10 9 M BH (Kormendy et al. 1996). 



This paper is organized as follows. In Section 2 we de- 
scribe our models. We then investigate the photometric ob- 
servables of the nuclear discs in Section 3, with attention to 
their effect on the surface brightness profiles and isophotes. 
In Section 4 we construct both the even and the odd part 
of the DF. We study the stability of the resulting models in 
Section 5. In Section 6 we discuss the kinematic signatures 
of nuclear discs and the influence of seeing convolution. We 
investigate the influence of a nuclear black hole in Section 
7, and summarize our results in Section 8. 

2 MULTI-COMPONENT MODELS 

We construct galaxy models consisting of a thick disc em- 
bedded in an spheroidal body: the 'bulge'. In addition, we 
allow a black hole to be present at the centre of the galaxy. 
The potential of the model can therefore be written as 

*(#, z) = ip bu i se (R, z) + tp disc (R, z) + iPbh(R, z). (2.1) 

We denote the total potential by ^(R, z), and use %j){R,z) 
to indicate the potential of a separate component. We use 
the relative potential Vl> = — <!> instead of the potential <3>. 
Therefore stars at rest in the centre of the potential well 
have the maximum energy. 

We define a disc to be nuclear if its horizontal scale 
length Rd is smaller than half the core radius Rb of the 
bulge. This is defined as the inner core- or break radius of 
the bulge component which ~ 2" — 3" at Virgo (Ferrarese 
ct al. 1994), and should not be confused with the 'the Vau- 
couleurs' effective radius which is much larger. Although 
this paper only focuses on such nuclear discs, the method 
discussed here is also applicable to larger discs, and even to 
SO galaxies, by simply increasing the ratio Rd/ Rb- 

Throughout this paper we take the mass-to-light ratio 
M/L to be equal to one, so that, e.g., the projected surface 
brightness is identical to the projected surface density. 

2.1 The bulge component 

We describe the spheroidal component by the so-called (a, 
(5) models, which have a density distribution given by 

^) = ^(g) a (l + ^)", (2.2) 
where 

m = \l R2+ w (2 ' 3) 

Here Rb is the 'break-radius' of the bulge, and q b is the 
flattening of the bulge. For a = the central density is 
finite and equal to po,f>- When a < the density profile 
has a central cusp. These models have proven to accurately 
fit the surface brightness distribution of elliptical galaxies 
(e.g., van der Marel et al. 1994, Qian et al. 1995 (hereafter 
QZMH), van den Bosch & van der Marel 1995). QZMH 
discuss these models in more detail. 

At large radii the density (2.2) falls off proportional to 
m Q+2/3 colder only models with a + 2/3 = —4, so that 
the projected surface density at large radii falls off as R~ 3 . 
In this case the total mass is finite, and given by 

M = 2^0,6^5(1 + 1,1), (2.4) 
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where B is the beta-function. These models are similar to 
the set of models discussed by Dehnen & Gerhard (1994), 
except that their density has (1 + m/Rt) 213 as second factor 
rather than (1 + m 2 / Rl) 13 . 

QZMH showed how the classical double quadrature ex- 
pression (e.g., Chandrasekhar 1969) for the potential of any 
axisymmetric system in which the density is stratified on 
similar concentric spheroids, i.e., p = p(m ), can always be 
rewritten as a single quadrature. For our particular choice of 
a + 2/3 = —4, however, the inner integration of the classical 
formula can be carried out explicitly, and therefore leads di- 
rectly to a simple single quadrature, which is evaluated more 
efficiently. According to the classical formula, the potential 
that corresponds to (2.2) is given by 



ip(R, z) = V>o,f> — 7rGqt 



where 



F(t) = / p(m' 2 )dm' 2 . 



F(t)di 



(T + l)^/r + q 2 



(2.5) 



(2.6) 



Substitution of the density (2.2) in the above integral, and 
changing to the integration variable x = ^—^577^^, yields 



F(t): 

with 

t = 



-p , b R 2 b ^ +2) (l+t 2 /R 2 f +1 



1+/ 



£2(0+1) 



R 2 

T + l 



(2.7) 



(2.8) 



This result is valid when (3 < — 1. Substitution of expression 
(2.7) in the integral (2.5) then gives the potential as a single 
quadrature. The central potential tp 0tb can be evaluated 
explicitly, and is given by 



tpo,b = 



~GM 



(l+/3)i? 6 B(f + §,!) yy—-^ 



(2.9) 



2.2 The disc component 

Recent work by, e.g., Scorza & Bender (1995), indicates that 
the small discs embedded in ellipticals form a smooth tran- 
sition from SOs to galaxies with smaller disc-to-bulge ra- 
tios. We therefore assume that our nuclear discs have a 
surface brightness distribution that is similar to the discs 
of SOs and spirals. Most studies have considered models 
with infinitesimally thin discs. Although such discs provide 
a good approximation for some purposes, they are not suffi- 
cient for our needs, because they are maximally cold in the 
perpendicular (z-) direction. Real discs have non-zero z- 
motions that contribute to the line-of-sight velocities when 
observed at any inclination other than edge-on. Further- 
more, infinitesimally thin models are unphysical in that the 
derivative dp/dz is infinite at z = 0. We therefore con- 
sider thick discs. Observations have shown that the surface 
brightness of discs is well approximated by an exponential 
profile, both in the radial and the vertical direction (e.g., 
Freeman 1970; Aoki et al. 1991). 



Few treatments of potential-density pairs of thick discs 
are available in the literature. Kuijken & Gilmore (1989) 
and subsequently Cuddeford (1993) presented a very sim- 
ple method of expanding infinitesimally thin disc potentials 
to a thick disc with arbitrary vertical density distribution 
with constant scale height. Consider an infinitesimally thin 
disc, whose surface density is p(R,z) — f(R)S(z), and with 
corresponding potential ip(R,z) — g(R,z). Here / and g 
are arbitrary functions. One can thicken this disc by replac- 
ing the Kronecker delta function S(z) by any other func- 
tion h(z) that describes the vertical density distribution. 
At each z 1 parallel to the equatorial plane the density dis- 
tribution can again be written with the delta function, i.e., 
p(R,z') = f(R)5(z — z'), and for this density distribution 
the corresponding potential is simply ip{R, z') = g(R,z — z'). 
The potential of the thickened disc can be built up of an infi- 
nite number of such infinitesimally thin discs, each of which 
has a corresponding potential, properly weighted by the ver- 
tical density distribution. The potential 



00 

-/ 



g(R, z — z')h(z')dz' 



(2.10) 



follows straightforwardly. 

Using the scheme outlined above, we could construct 
the potential density pair of the double exponential disc, i.e., 
p(R,z) = Eo exp(— R/Rd) exp(— |a;|/^o)- The potential of 
the infinitesimally thin case is already a single quadrature, so 
that equation (2.10) becomes a double integration. For the 
special case of a double-exponential disc this can be reduced 
to a single quadrature (Kuijken & Gilmore 1989). This is 
an important advantage, since it speeds up the numerical 
calculations considerably. However, the double-exponential 
has the unphysical characteristic that the derivative dp/dz 
is discontinuous at z = 0. Furthermore, the slow oscilla- 
tory behavior of the Bessel functions that occur in the inte- 
grand makes the numerical evaluation of the potential and 
its derivatives tedious. 

We therefore decided to use an 'exponential spheroid', 
first introduced by Kent, Dame & Fazio (1991), to model 
the disc. This model, like the double-exponential disc, has 
an exponential surface brightness along both the radial and 
the vertical direction. 

Consider a sphere whose projected surface brightness is 



E(i?) = E exp{-R/R d ) 



(2.11) 



The corresponding density distribution follows upon solving 
the Abel integral equation 



E(i?) 



OO 



prdr 



vV 2 - R 2 : 



(2.12) 



which gives 

p(r)=p 0ld K {r/R d ), (2.13) 

where po,d = ^o/irRd, and Ko is a modified Bessel function. 

We replace the sphere by a sp heroid of fl attening qd, 
and introduce the spheroids m = R 2 + z 2 Jq%. This leads 
to axisymmetric models that project to a surface brightness 
that is exponential along any axis, i.e., 



= ^S exp(- J R7i? d ). 



(2.14) 
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Here R* = ^x 2 + y 2 /q2 and q d is the projected axis ratio, 
which is related to the intrinsic flattening and the inclination 
angle i through 

q'i = cos 2 i + q 2 d sin 2 i. (2-15) 

When qd <C 1 these models represent thick exponential discs. 

The potential of the above discs can be evaluated by 
means of equation (2.5). Substitution of equation (2.13) in 
equation (2.6) gives 

F(t) = -2po,di?d^i(^-), (2.16) 
so that the potential of the exponential spheroids is given 

by 



MR,*) 



GM D 



oo 

/ 



ttfl(£)dT 



(2.17) 



but with q b replaced 



where t is again given by equation (2.1 
by q d - 

The potential can be written as a single quadrature with 
the modified Bessel function K\ in the integrand. Since 
these functions decay exponentially and are positive defi- 
nite, their numerical evaluation is much simpler than in the 
case of the double-exponential disc mentioned above. Fur- 
thermore, the derivative dp/dz is continuous at z = 0. 

The central density of the exponential spheroid is infi- 
nite, resulting from the fact that Yaa x ^oKo(x) = —\n(x). 
However, the total mass is finite; M D — 2-K 2 q d po, dRd- The 
central potential is given by 



i>0,d 



2GM D arcsin(yr 



•q3) 



TlRd 



(2-18) 



2.3 The parameters 

In order to constrain the large parameter space of our multi- 
component models we fix the values of a number of the pa- 
rameters. Since we are mainly interested in the central parts 
of the models the only bulge-parameter that we vary is the 
cusp steepness a. The other parameters (R b , q b , and p ,&) 
are kept fixed. The value of /3 is set by our requirement that 
a + 2/3 = — 4. We choose a flattening of 0.7 for the bulge. 

For the disc we choose a fixed flattening of qd = 0.1. 
This is consistent with the observed flattening of several 
nuclear discs (van den Bosch et al. 1994). Instead of param- 
eterizing the disc by its mass Md and scalelength Rd, we use 
the ratio Rd/Rb and the disc-to-bulge ratio A. We restrict 
ourselves to models with R d /Rb = 0.2. This value is chosen 
to match approximately the values observed in a number of 
Virgo ellipticals (Scorza & van den Bosch, in preparation). 
Furthermore, taking R d /R b — 0.2 ensures that our choice to 
restrict ourselves to models with a + 2/3 — —4 does not influ- 
ence the results. This is due to the fact that at radii where 
the slope of the density distribution of the bulge starts to 
depend on the particular choice of a + 2/3 (i.e., at R > Rb), 
the contribution of the nuclear disk is already negligible. 

We define A as the ratio of the total luminosity of the 
disc to the total luminosity of the bulge component. We 
assume that the mass-to-light ratios of both components are 
equal, and then have 



Table 

Model 

1.0 

1.1 

1.2 

1.3 

1.4 

2.0 

2.1 

3.0 

3.1 

3.2 

3.3 

3.4 

3.5 

3.6 

4.0 

4.1 



Parameters of 

a A 
0.0 
0.0146 
0.0146 
0.0146 
0.0146 

0.0 
0.0207 

0.0 
0.0303 

0.0 
0.0303 

0.0 
0.0303 
0.0303 

0.0 
0.0444 



the models 



0.0 
0.0 
0.0 
0.0 
0.0 
-0.5 
-0.5 
-1.0 
-1.0 
-1.0 
-1.0 
-1.0 
-1.0 
-1.0 
-1.5 
-1.5 



A 
0.0 
1.0 
1.0 
1.0 
1.0 
0.0 
1.0 
0.0 
1.0 
0.0 
1.0 
0.0 
1.0 
1.0 
0.0 
1.0 



M BH /M bul , 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 
0.00303 
0.00303 
0.0303 
0.0303 

0.0 

0.0 

0.0 



fo(E,L z 
ST 
ST 
NR. 
NR. 
CR 
ST 
ST 
ST 
ST 
ST 
ST 
ST 
ST 
CR 
ST 
ST 



+10 
+10 
+10 
-10 
+10 
+10 
+10 
+10 
+10 
+10 
+10 
+10 
+10 
+10 
+10 
+10 



M, 



disc 



Mbulgo 



Po,d /Rd\3 



B(f + §,i)©Po,* v JV 



(2.19) 



For each model, the cusp-steepness a, the disc-to-bulge ratio A, 
the ratio of surface brightness for disc and bulge at R = R d (edge- 
on) A, and the ratio of BH mass over bulge mass are shown. In 
addition, the method and the value of the parameter a used to 
calculate the odd part of the DF are given. Here ST stands for the 
'standard' models, NR for 'non-rotating bulge' models, and CR for 
'counter-rotating' models (see Section 4.2). 

Figure 1 shows a contour-plot of the potential and density 
in the meridional plane of a model with a = —1.0, and 
A = 0.0303. Although the presence of the nuclear disc is 
hardly visible from the potential contours, the flattened den- 
sity contours in the centre clearly reveal its presence. 

We also define the ratio of the surface brightness of disc 
and bulge at R = R d , for an inclination angle of i = 90°, 

A _ 90 (R d ,0) 

We choose the disc-to-bulge ratio A of our models in such 
a way that A = 1.0. In this way, the relative contribution 
of disc and bulge to the velocity profiles at R = Rd is equal 
(for i = 90°). 

Table 1 summarizes the parameters of the entire set of 
models. 



3 PHOTOMETRIC OBSERVABLES OF 
NUCLEAR DISCS 

We first investigate the photometric characteristics of nu- 
clear discs embedded in a spheroidal body. We project 
the models defined in Section 2 at different inclination an- 
gles i, and subsequently seeing-convolve them with a Gaus- 
sian Point Spread Function (PSF). The parameter x = 
FWHM/7? d gives the FWHM of the PSF expressed in units 
of the horizontal disc scale length Rd- We scale the galaxy 
such that the FWHM corresponds to 1" in all cases, and 
take our pixels to be 0.2" x 0.2". After seeing-convolution 
we determine, as a function of radius, the ellipticity and 
the amplitude of the cos(4#)-term in the Fourier expansion 
that describes the deviation of the isophote from a pure el- 
lipse (e.g., Lauer 1985, Jedrzejewski 1987). If this amplitude 
is positive it indicates that the isophote is 'discy', whereas 
negative amplitudes imply 'boxy' isophotes. 

Figure 2 presents the results. Since we define the 
FWHM of the PSF to equal 1", changing \ corresponds 
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to changes in the scale of the galaxy, i.e., its distance. We 
plot all parameters out to 10". The presence of the nuclear 
disc is evident from the central increase in ellipticity, and the 
disciness of the isophotes in the central region. When the 
bulge harbours a steep cusp the main photometric signature 
of the nuclear disc is the disciness of the isophotes, as can be 
seen in the panels on the left. The panels in the middle of 
Figure 2 show that the nuclear discs become photometrically 
undetectable if the FWHM of the PSF exceeds 2-3 times 
the disc scale length. Furthermore, as was already shown 
by Rix & White (1990), the discs become photometrically 
undetectable for i < 60° (panels on the right). 

When trying to decompose the surface brightness distri- 
bution in a disc and a bulge component from the photometry 
alone, a number of indicators can be used. The radius at 
which the cos4# profile reaches its maximum is strongly cor- 
related with the scale length Rd of the disc (Scorza &: Bender 
1995). The left panels in Figure 2 illustrate that the central 
ellipticity profile contains information on the cusp steepness 
q. The amplitude of the central changes in the ellipticity and 
cos4# profiles are related to the disc-to-bulge ratio (or equiv- 
alently, the central surface brightness of the nuclear disc). 
All these indicators also depend on the inclination angle i. 
Scorza & Bender (1995) showed that the cos6# parameter 
can be used to constrain i, under the assumption that the 
disc is infinitesimally thin. However, there is a degener- 
acy between disc thickness and inclination angle, because 
our discs are stratified on similar concentric spheroids. This 
will lead to some non-uniqueness for the disc and bulge pa- 
rameters, but the solution space will be small for sufficiently 
inclined galaxies. 

4 CONSTRUCTION OF THE F(E,L Z ) MODELS 

Any collisionless system can be fully specified by its distribu- 
tion function, which, as stated by the Strong Jeans Theorem, 
depends on the phase-space coordinates only through its iso- 
lating integrals of motion. Most stellar orbits in an axisym- 
metric potential admit three such isolating integrals, namely 
the energy E, the angular momentum along the rotation axis 
L z , and some third integral 73. Except for Stackel poten- 
tials (e.g., de Zeeuw 1985), the non-classical third integral is 
generally not known analytically, and may not exist for all 
orbits. This is the reason that most axisymmetric models 
constructed to date have distribution functions that depend 
only on the two classical integrals of motion; / = f(E,L z ). 
Although such models often predict too much motion on the 
major axis due to the fact that the velocity dispersions in 
the R- and ^-directions are everywhere equal, or = o z , they 
provide reasonable approximations to the dynamics of some 
galaxies (van der Marel 1991; QZMH; Dehnen 1995). here 
we also limit ourselves to two-integral models. We discuss 
the limitations of this approach in Section 8. 

4.1 The even part of the distribution function 

The density that corresponds to the combined potential 
ty{R, z) of bulge, nuclear disc, and BH, depends only on 
the part of the DF that is even in L z : 

y R A /2(*-E) 

f e (E,L z )dL z . (4.1) 



P=^lAE 



R 



For L z = the inverse of this equation is given analytically 
by Eddington's formula 



fe(E,0) = 



d 2 p d* 



d 2 * y/E- * 



Je W*/ 



(4.2) 



but the evaluation of f e {E, L z ) for other values of L z from a 
given p(R, z) is not easy. Lynden-Bell (1962), Hunter (1977), 
and Dejonghe (1986) used different transformation methods 
to evaluate f e (E,L z ) for a number of models. However, 
since all these methods require the analytical knowledge of 
p(^,R), they are not widely applicable. 

Recently more general schemes have become available 
for the evaluation of f e (E,Lz)- Hunter & Qian (1993, 
hereafter HQ) found a contour integral expression, which is 
the generalization of Eddington's formula. Other methods 
for generating f e (E,L z ) were developed by Dehnen (1995), 
Magorrian (1995), and Kuijken (1995). None of these meth- 
ods requires that the density can explicitly be written as 
a function of so they allow the evaluation of f e (E,L z ) 
for any p(R,z). QZMH and also Dehnen (1995) success- 
fully used these methods to construct a two-integral model 
of M32 with a central black hole that provides an accurate 
fit to the observed velocity profiles. 

We use the HQ-method. The HQ solution for the even 
part f e (E,L z ) of the DF is 



fe(E,L z ) = 



[*e„„(£)+] 



V8n 2 



I 



Pn 



2(C - E) 



(4.3) 



Here the integral is along a complex contour in the plane of 
the complex potential £. The two subscripts of pn denote 
the second partial derivative with respect to the first argu- 
ment (£). The tilde denotes the complex continuation of the 
function pn, which is given by 



P22(R 2 ,z 2 ) p 2 (R 2 ,z 2 )^ 22 ((R 2 ,z 2 



[<K 2 (R 2 



[* 2 (i? 2 



(4.4) 



Each subscript 2 denotes the partial derivative with respect 
to z 2 . The HQ paper gives further details on the method. 
The complete procedure of how to calculate the complex 
contour integral is discussed in detail in QZHM. 

We split f e (E, L z ) in two parts; a bulge part fe(E, L z ), 
and a disc part f^(E,L z ). Each of these two terms of 
f e (E,L z ) is calculated by solving equation (4.3), replacing 
P = Pbuige + Pdisc with pbuige and p disc respectively. The 
potential remains the total potential * = i/Wge + V>disc- 

We evaluate f e (E, L z ) = f e {E, L z )+ff{E, L z ) on a 40 x 
40 grid in the region of the (E, L 2 )-plane that corresponds 
to bound orbits in the potential VF This region is bounded 
by L z > 0, E < *oo, and by the curve L z < Z/ Z , max (-E). 
Here L Ztln ^(E) is the maximum allowed value of L z for 
energy E, and corresponds to the angular momentum of a 
star with energy Bona circular orbit. For the 40 grid points 
that have zero angular momentum we compare the result 
of the complex contour integration with the value given by 
Eddington's formula (4.2), to check that the integral on the 
chosen contour converges, and whether any adjustments of 
the contour are necessary. We check for self-consistency by 
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calculating p(R,z) from the DF, using equation (4.1). Bi- 
cubic spline interpolation is used to interpolate between grid 
points. The models are self-consistent to a high degree of 
accuracy (~ 10 ). 

Figure 3 shows log 10 [/ e (S, L z )] as function of E/^o and 
rj 2 = [L z /L Zz max (E)] 2 for four different models: model 1.0 
(upper left), model 1.1 (lower left), model 3.0 (upper right), 
and model 3.1 (lower right). See Table 1 for the parameters 
of each of these models. The 40 x 40 grid is also shown. In 
order to properly sample the DF, the grid-density increases 
with energy and angular momentum. The presence of a 
central cusp in the bulge density is evident from a strong 
rise in the centre (i.e., at large E/^q). The presence of the 
nuclear disc also induces such a strong rise. Furthermore, 
close to the centre the DF is peaked towards high angular 
momentum, indicative of the rotational support of the disc. 

4.2 The odd part of the distribution function 

Whereas the even part of the DF generates the density of 
the system through equation (4.1), the odd part of the DF 
specifies the mean streaming motion: 



47T 

R 2 



g, Ry/2(V-E) 

dE J f a {E,L z )L z AL z 
o o 



(4.5) 



Therefore, only if in addition to p(R,z), the entire mean 
streaming of the stars is known a priori, can one obtain 
fo(E,L z ) by inversion of equation (4.5). Since we are pri- 
marily interested in studying the observable properties of 
nuclear discs for different amounts of rotational support, we 
exploit the freedom to specify fo(E, L z ) under the constraint 
that f(E,L z ) ee f e (E,L z )+f (E,L z ) be positive. 

We consider three different approaches for specifying 
f e (E, L z ). In the first approach we define a function h a (rf), 
where n = L z /L Zjmax (i5), and simply take 



U{E,L z )=h a (ri)U{E,L z ), 



(4.6) 



By doing so we consider the entire system as a one- 
component model, i.e., although the potential is provided by 
the disc and the bulge, the two components are dynamically 
coupled to each other since we do not treat the streaming 
motions in them separately. In order to limit the number of 
free parameters we consider a f (E,L z ) that is fully spec- 
ified by only one parameter. Dejonghe (1986) showed that 
there is such a functional form which has the advantage 
that it maximizes the entropy of a two-integral axisymmet- 
ric system. We follow van der Marel et al. (1994), and use a 
modified version of this parameterization, which is given by 

r tanh(a?7/2) / tanh(o/2) (a > 0) 

ha{v) = \v (a = 0). (4.7) 

[ (2/a)arctanh[?7tanh(a/2)] (a < 0) 

The other advantage of this parameterization is that these 
models will automatically be physical (i.e., have positive 
DF) as long as f e (E,L z ) > 0. The free parameter a deter- 
mines the amount of rotation, whereby the model of max- 
imum rotation has a = co. Van der Marel et al. (1994) 
and also QZMH used an extended version of this param- 
eterization in which they included another parameter that 
determines the fraction of stars on clockwise, circular orbits 



in the equatorial plane. Although this allows the study of 
more general models, it invokes yet another free parameter. 
Therefore, we restrict ourselves to models in which this frac- 
tion is unity, i.e., in which all the stars on circular orbits in 
the equatorial plane have the same rotation direction. 

The second approach allows us to construct models 
with a non-rotating bulge and a rotating disc. We calcu- 
late fe(E, L z ) and f^(E, L z ) separately, and then define the 
odd part of the distribution function of disc and bulge sepa- 
rate from each other. This gives models that consist of two 
(dynamically) decoupled components. We consider cases in 
which fa(E,L z ) — (no rotation of the bulge component), 
and 



ft(E,L z )=ha(ji)ft(E,L z ). 



(4.8) 



Here h a (rj) is defined by equation (4.7). 

Finally, in our third approach for specifying f (E,L z ) 
we construct models with a counter-rotating disc by simply 
taking 



fo(E,L z ) 
and 

fo(E,L z ) 



+h a ( V )f°(E,L z ), 



-ha(v) fe(E,L z 



(4.9) 



(4.10) 



Throughout the remainder of this paper we will refer to mod- 
els constructed with the first method (i.e., with dynamically 
coupled components) as 'standard models'. The models for 
which the odd part of the distribution function of the bulge 
is zero as 'models with non-rotating bulge', and the models 
whose odd part results in counter-rotating discs as 'counter- 
rotating models'. 

4.3 The region of physical models 

QZMH showed that two-integral oblate (a, /3)-models with a 
nuclear BH are unphysical (i.e., have a region in phase-space 
where the DF is negative) when a > —0.5. Therefore, it is to 
be expected that models with very compact, massive nuclear 
discs embedded in an (a, /3)-bulge, with a > —0.5 will also 
be unphysical. Since a nuclear disc is a cold component (see 
also Section 5.2) located at the centre of the galaxy, most 
stars in the centre will be on nearly circular orbits. In the 
derived DF this translates into small weights for the more 
radial orbits close to the centre. This can clearly be seen 
from the DF of model 1.1 in Figure 3. The DF is peaked 
towards large energies (i.e., the centre of the potential well) 
and high angular momentum. A local minimum close to the 
centre at L z = is also visible. 

The value of the DF at this local minimum depends 
on the disc-to-bulge ratio of the system. The models will 
become unphysical if A = Mdisc/Afb u i g c is too large, since 
this will bring about a region in phase-space were f e (E,L z ) 
is negative. If A is increased further the nuclear disc be- 
comes completely dominant over the bulge. In those cases 
the DF is positive again over the entire phase space. So 
as function of Rd there is an interval of disc-to-bulge ratios 
for which the DF has a local minimum with negative val- 
ues. Since the angular momentum at this local minimum 
is zero, one can test whether a certain (Rd/Rt, A)-model 
is physical by examining whether f e (E,L z = 0) is positive 
for all E < V&o- In order to determine this interval we pro- 
ceed as follows. For each model (Rd/Rb,A) we numerically 
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evaluate E CI it for which f e (E,L z — 0) is (locally) minimal, 
using Brent's method for the minimization. Subsequently, 
we check whether f e (E crlt ,L z = 0) < 0. We repeat this pro- 
cedure in a bisective way until the borders of the interval are 
known to a certain accuracy. The values of f e (E,L z — 0) 
are evaluated using Eddington's formula (4.2). 

The results for the case with a — and a — —0.25 
are shown in Figure 4, where we plot log 10 (A) as function 
of \og 10 (R d /R b ). The hatched areas indicate the regions 
where the models are unphysical. Only very small, massive 
nuclear discs give rise to unphysical two-integral DFs. As 
expected, the region of unphysical models becomes smaller 
when the cusp steepness increases (i.e., when the value of a 
decreases) . 



5 THE STABILITY OF THE NUCLEAR DISCS 

Our f(E, L z )-models are of little practical value if they are 
unstable, and hence we study their stability in this section. 



5.1 Local stability 

The best known stability criterion for discs is the Toomre 
(1964) criterion. This criterion ensures that an infinitesi- 
mally thin disc is locally stable against axisymmetric per- 
turbations as long as 



Q 



7 GE 



> 1. 



(5.1) 



Here or is the radial velocity dispersion, E is the surface 
density of the disc, 7 is a constant that equals 3.36, G is the 
gravitational constant, and k is the epicyclic frequency. 

In an axisymmetric system with distribution function 
f e (E,L z ) we have or — a z = a everywhere. From the 
Jeans equations it follows that (e.g., Hunter 1977) 



a\R,z)= ' 



The surface density of our disc is given by 

00 

T,(R) = / p disc (R,z)dz. 



(5.2) 



(5.3) 



A problem that arises is that we are considering a thick disc. 
Only a limited amount of work has been done on the stability 
of thick discs, but the little work that has been done (e.g., 
Shu 1968; Vandervoort 1970) indicates that disc thickness 
has a stabilizing effect. This means that the factor 7 = 3.36 
derived by Toomre for the infinitesimally thin case decreases 
for thick discs. However, it is not known what the value of 7 
is for an arbitrary vertical density profile. Therefore, instead 
of presenting Q(R), we show the radial profile of yQ(R). 

The results are shown in Figure 5, where we give 7Q 
as function of radius for four models that differ only in cusp 
steepness. The hatched region indicates the regime where 
infinitesimally thin discs are unstable against axisymmetric 
perturbations (i.e., 7Q < 3.36). In the case of a thick disc, 
the upper limit of this regime will move slightly downwards. 
The presence of a cusp stabilizes the disc in the inner region. 
Only the nuclear disc in the model without cusp (a = 0.0) 



is mildly unstable in its inner region. Figure 6 shows the 
region in the (a, A)-plane of models that are stable over the 
entire radial range (hatched region) , where we have assumed 
that stability is achieved for 7Q > 3.0. Although this value 
is uncertain, as we have seen, a 10% change will cause only 
a very small shift in the location of the boundary curve. It 
is again evident that a cusp has a strong stabilizing effect 
on a nuclear disc. 

Although the surface density of the disc falls off 
strongly, the velocity dispersion decreases slowly, due to 
the fact that there is a large amount of matter (from the 
bulge-component) beyond several scale lengths of the disc. 
Therefore, the outer parts of the disc are stable in all cases. 

Figure 7 shows 7Q as function of radius for the cusp-free 
model with A = 1.0 (model 0.1), for three different fiatten- 
ings of the nuclear disc. Thickening the disc, i.e., increasing 
qd, has a stabilizing effect. This can also be understood 
intuitively: the thickness of a disc is closely related to the 
amount of random motion of the stars in the ^-direction. 
Since our models have DFs that depend on only two inte- 
grals of motion, we have or — a z . Therefore, thickening the 
disc will increase the radial velocity dispersion in the disc, 
and therefore the parameter Q. The dependence of k on qd 
is weak. 



5.2 Global stability 

Besides local instabilities, discs can also suffer from global 
instabilities. The best known is the bar-instability, which 
has proven to play an important role. Ostriker & Peebles 
(1973) calculated that a disc is globally stable to bar- like 
modes when the total kinetic energy of rotation T is less 
than ~ 14% of the gravitational energy \W\. However, 
this criterion is merely a good rule of thumb, rather than a 
proven physical theorem. Athanassoula & Sellwood (1986) 
showed that, even more effectively than a halo component, 
random motions in the disc can stabilize it against bar form- 
ing modes. They showed that as long as the mass-weighted 
value of Q averaged over the central region exceeds a value 
of ~ 2 — 2.5, the entire disc is stable. This is the case for 
all of the models listed in Table 1. 

Toomre (1981) argued that the origin of the bar form- 
ing instability is feedback to the swing amplifier via waves 
that pass through the centre. The occurrence of swing am- 
plification depends on the value of 



X 



n 2 R 



2m^' ^ 
where m = 2 in the case of bi-symmetric bars. Toomre 
showed that as long as Q < 3 and X < 3 the gain of the 
swing-amplifier is sufficient to form bars. We find that X 
increases strongly with radius, but X < 3 for small R. Al- 
though this implies that swing amplification can occur in 
the very centre of the nuclear discs, one can ensure stability 
against bar formation by cutting off the path to the centre 
by means of an inner Lindblad resonance. The models with 
a cusp (a < 0) have Q — — > +00 for R — » 0, and will 
therefore have an ILR for any pattern speed of the bar. An 
ILR is absent only when a — and high pattern speeds. In 
that case the swing amplified feed-back loop is open in the 
centre (X < 3), but as there is a large amount of random 
motion in the centre (i.e., Q > 3) the growth rate for a bar 
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mode will be insignificant. Therefore, we expect that the 
nuclear discs in all the models listed in Table 1 are stable 
against bar formation. 

The above results are valid for a bulge with axis ratio 
of 0.7. N-body experiments carried out by Dehnen (priv. 
comm.) indicate that non-rotating and rotating (a, j3)- 
models with axis ratios of 0.5 or smaller are unstable to 
bar formation. 



5.3 Continuation to larger discs 

We have restricted ourselves to nuclear discs with Rd/Rb = 
0.2, but the method outlined above is also applicable to 
larger discs. We have investigated the stability of our disc- 
models for a large range in Rd/Rb- The general tendency 
is for larger discs to be more stable. We have calculated 
the maximum disc-to-bulge ratio for which, at given cusp- 
steepness a and disc scale length R d , ^yQ(R) > 3.0. From 
equation (2.19) it follows that Eo oc A/R^. Therefore, for 
given disc-to-bulge ratio and disc-scale length we can, using 
an arbitrary scaling, calculate the central surface brightness 
(M) of the disc (in magnitudes per arcsec ) . 

Figure 8 shows the curve of fxo as a function of 
\og 10 [Rd/Rb] for discs with maximum disc-to-bulge ratio so 
that the entire disc has 7Q > 3. These curves depend on 
the cusp steepness a, but only for small discs. This Figure is 
remarkably similar to the relation between /io and Rd of ob- 
served discs (i.e., compare to Figure 17 of Scorza & Bender 
1995). This might indicate that discs build up their mass 
until they become marginally stable. 



6 THE KINEMATICS OF NUCLEAR DISCS 

Our models are fully specified by the four parameters a, 
A, a, and Mbh- In addition, we must choose the method 
used to calculate the odd part of the DF. In this section we 
restrict ourselves to cases with Mbh = 0. The influence of 
a nuclear BH in the centre of these models is the topic of 
Section 7. 



6.1 Calculation of the velocity profiles 

Once we have determined the complete phase-space density 
of our models, we can calculate every observable property. 
In particular, the entire VP can be derived by integrating 
the DF along the line of sight: 

VP{v t ,;x',y') = ±J J J f(E, L z )dv x ,dv y ,dz' . (6.1) 

-E>*oo 

Here we adopt the Cartesian coordinate system of an ob- 
server (x , y' , z'), where the «'-axis lies along the line of sight, 
and x' is along the apparent major axis of the galaxy pro- 
jected on the sky. 

Upon defining a polar coordinate system (v± , if) in the 
(v x i,v y i) plane, where v x i = v±cos(p, and v y i = wxsiny, 
equation (6.1) can be written as 



22 z 1 ir 

VP(v z ,;x',y') = ± J dz' J dvl J f(E,L,)d<p. (6.2) 



The energy is given by 

E = *(x',y\z')-±(v z ,+vl), 

and the angular momentum follows from 



+vi_cos(pyJ (— y'cosi + 2'sini) 2 + a;' 2 cos 2 i 



(6.3) 



(6.4) 



For any given v z i the boundaries zi and 22 are determined 
by solving ^f(x',y',z') - \v\, = 0. 

For each model we calculate the VPs at a number of 
(a/, j/')-points. In all cases discussed here we take an in- 
clination angle of 90°, i.e., we assume edge-on observation. 
Figures 9a and 9b compare the VPs of model 1 (a. — 0.0) 
at seven positions along the major axis (x' = R, y' = 0), us- 
ing four different odd parts of the distribution function: the 
first (left panels in Figure 9a) has f (E,L z ) specified by the 
'standard' method (Section 4.2). The panels on the right in 
the same figure correspond to a 'non-rotating bulge' model 
with a — +10, whereas taking a = —10 results in the VPs 
shown in the left panels of Figure 9b. Finally, in the right 
panels of Figure 9b, the VPs are shown that correspond to 
the case of a counter-rotating disc, i.e., where we have re- 
versed the sign of the odd part of the DF of the disc with 
respect to that of the bulge. 

Even though the dynamics of bulge and disc are coupled 
in the 'standard' models, the VPs nevertheless show two 
very distinct components. The solid triangles in Figure 9 
indicate the circular velocity of the models. As can be seen, 
the disc rotates with almost this circular velocity. This is 
true for both the 'standard' models, the 'non-rotating bulge' 
models, and the 'counter-rotating' models. 

As discussed in Section 4.2, we have constrained the 
fraction of stars on clockwise circular orbits in the equatorial 
plane to be unity. However, not all stars in the disc are on 
such circular orbits. Therefore, upon taking a <C in the 
odd part of the DF, a considerable fraction of stars in the 
disc will be on anti-clockwise orbits that are not (perfectly) 
circular. This results in double-peaked VPs, as can be seen 
in the left panels of Figure 9b. 



6.2 Velocity profile analysis 

In order to allow for an easy comparison between different 
models, and to be able to extract physically meaningful pa- 
rameters, we need to find a way of parameterizing the VPs. 
It has become customary to expand the VPs in a Gauss- 
Hermite series (van der Marel & Franx 1993; Gerhard 1993). 
Although this approach has proven useful when parameter- 
izing VPs that deviate slightly from a Gaussian, it is in- 
appropriate for describing the VPs presented above, since 
they strongly deviate from Gaussians. After seeing convo- 
lution with a FWHM that exceeds a few disc scale-lengths 
the Gauss-Hermite parameterization can be used, since the 
convolution results in sufficiently smooth VPs. However, 
for smaller seeing FWHM the convolved VPs still deviate 
strongly from a Gaussian. 

We therefore opted to use the real moments to quantify 
the VPs. The n th moment of the normalized velocity profile 
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VPo(w) is defined as 

oc 

J VPo(w) / dv. (6.5) 

— oo 

The mean streaming rotation velocity is simply the first mo- 
ment V ro t = /J-i, the velocity dispersion a = y/ (12 — n\ , 
the skewness S — (/is — 3/tt2/ii + 2/il)/a 3 , and the kurto- 
sis K = (/t*4 — 4fi 3 fj,i + 6fi2fJ,i — 3fif)/a 4 . 

Figure 10 shows these parameters for models 1.1 - 1.4 
whose VPs are given in Figures 9a and 9b. The presence 
of the nuclear discs brings about VPs that deviate strongly 
from a Gaussian. This is evident from the large values of 
both S and JC, which are zero and three respectively for a 
Gaussian. The disc dominates the kinematics in the centre, 
and causes a strong drop of the velocity dispersion inwards 
of ~ 3i?d- The rotation curve of model 1.4 clearly reveals a 
counter-rotating core. This is also evident from the change 
of sign of the skewness S around ~ 2.5Rd- 

The real moments of the VPs describe the dynamics of 
the entire system. If one is interested in the dynamics of 
disc and bulge separately, another parameterization of the 
VPs is required. Given the clear two-component character 
of the VPs, it is convenient to use the Levenberg-Marquardt 
method to fit a double Gaussian to the derived VPs. Under 
the assumption that the two separate components of the VPs 
correspond to the disc and bulge components of the model, 
and assuming that both are Gaussian, the double-Gaussian 
fit measures the kinematics of disc and bulge separately (Rix 
& White 1992). There is no a priori reason why the VPs of 
both components will be Gaussian. However, increasing the 
number of parameters of the fitting function will result in 
a large degeneracy of the parameters. Furthermore, Scorza 
& Bender (1995) have shown that including deviations from 
a Gaussian only mildly influences the derived kinematics of 
bulge and disc. 

6.3 Observations of nuclear discs 

We now discuss the influence of seeing convolution on the 
observable properties of nuclear discs. We mimic long-slit 
observations where we take the slit to be aligned along the 
major axis, the slit-width to be equal to the FWHM of the 
seeing, and the pixels to be one-third of the seeing FWHM. 
We use the parameter \ to define the PSF FWHM, in units 
of the scale-length Rd of the disc, and subsequently scale 
the models so that the FWHM equals 1". Therefore vary- 
ing x is similar to varying the distance of the galaxy from 
the observer. We start by constructing a cube that con- 
tains the velocity profiles VP(iv; x' , y') calculated on a two- 
dimensional logarithmic grid on the sky (x' , y'). We then 
convolve the velocity profiles for each v z i with a Gaussian 
PSF and take the effects of slit-width and pixel size into ac- 
count. We use bi-cubic spline interpolations to interpolate 
between (x',y') grid points at a given velocity v z i. We use 
the approach taken by QZMH, who have shown that in the 
case of a Gaussian PSF and rectangular pixels, the entire 
seeing convolution and pixel averaging can be written as a 
double quadrature. After the convolution we parameterize 
the VPs. Since the convolved VPs deviate strongly from a 
Gaussian for small values of x, we use the moments defined 
in the previous section to describe the VPs. 



Figure 11 shows the results for four different values of 
the cusp steepness a. They include the case without nuclear 
disc (solid lines) as well as the case with nuclear disc with 
disc-to-bulge ratio such that A = 1.0 (dashed lines). The 
odd part is determined in the 'standard' way with a — 10. 
In all cases we take an inclination angle of i = 90°, and con- 
volve the model with a PSF with x — 1- The main effect 
of changing a is that the rotation curve of the entire galaxy 
becomes steeper, due to larger central mass. The kinematic 
signature of a nuclear disc is therefore most pronounced in 
models with steep cusps. The seeing convolution results in 
velocity dispersions that are almost equal for the cases with 
and without nuclear disc. The strong central decrease of 
dispersion (as can be seen in Figure 10) is no longer observ- 
able for a seeing with FWHM equal to the scale-length of 
the nuclear disc. The kinematic signatures of the nuclear 
disc are confined to higher rotation velocities and stronger 
deviations of the VPs from a Gaussian. 

Figure 12 shows in more detail the effect of the seeing 
convolution. Here we plot the moments of the convolved 
VPs for the model with a — —1.0, again both with and 
without nuclear disc (models 3.0 and 3.1 respectively), for 
four different values of x- Upon increasing x> the VPs of 
the models with and without nuclear discs become more and 
more similar. For x ~ 3 the only kinematic signatures of the 
nuclear disc that remain are a slightly larger maximum ro- 
tation velocity and central velocity dispersion. However, the 
differences are so small that it is unlikely that the presence 
of a nuclear disc can be inferred from kinematic observa- 
tions with a seeing FWHM exceeding ~ 3 times the disc 
scale-length. This is also the maximum x beyond which the 
disc becomes photometrically undetectable. Figure 13 sum- 
marizes these results for the velocity dispersion. We give the 
ratio of central velocity dispersion for the case with nuclear 
disc (A = 1) to the case without nuclear disc (A = 0) as func- 
tion of x f° r f° ur different values of cusp steepness. Increas- 
ing a results in an increase of the central mass, and therefore 
in an increase of the central rotation velocity. Increasing x 
results in an increase of the ratio <ro(A = 1)/<to(A = 0) due 
to seeing convolution of the rotation curve. 

A number of authors have suggested that the presence 
of a nuclear disc can mimic the presence of a BH, in that 
seeing convolution of the central rotation curve might re- 
sult in an increase of the central velocity dispersion. If the 
PSF FWHM is so large that the disc is photometrically un- 
detectable, the large central velocity dispersion could erro- 
neously be interpreted as due to a BH. However, we are 
unable to achieve a ratio of oo(A = l)/cr (A = 0) larger 
than 1.1. This is due to the fact that the PSF is circular 
compared to the highly flattened disc structure. Upon re- 
ducing the slit-width we only find a very mild increase of 
(To (A = l)/i7o(A = 0). The only way of achieving a large 
enough ratio so that the nuclear disc might mimic a BH 
is by strongly increasing a and the disc-to-bulge ratio A. 
However, stability arguments prohibit too large values of 
A. Furthermore, increasing A will make the nuclear disc 
detectable in the photometry. 

The entire analysis above is based on edge-on projec- 
tions of the models. Changing the inclination angle has 
a number of effects. First of all, since the disc is highly 
flattened, the observed rotation of the disc roughly scales 
with sini. Furthermore, decreasing the inclination angle 
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decreases the surface brightness of the disc relatively more 
than that of the bulge component. As a consequence, de- 
creasing i results in VPs that less and less clearly show a 
rapidly rotating component. We projected model 3.1 at a 
number of different inclination angles and found the kine- 
matic signatures of the nuclear disc to disappear for i < 60°. 

6.4 Counter-rotating discs 

In Section 4.2 we discussed how to define the odd part of 
the DF in order to achieve a model in which the nuclear disc 
counter-rotates with respect to the bulge. Several galaxies 
have been observed to harbour counter-rotating cores. Re- 
cent imaging with the Hubble Space Telescope has shown 
that in a number of cases there are indications for the pres- 
ence of a disc component in the centre (Forbes, Franx, & 
Illingworth 1995, Carollo et al. 1996). Therefore, counter- 
rotating discs seem a natural explanation for the observed 
kinematics. Here we investigate the kinematic signatures of 
such counter-rotating nuclear discs. 

We find that for x ~ 2 the VPs are double-peaked, 
clearly revealing the counter-rotation. In these cases one 
observes also a strong central decrease of the velocity dis- 
persion. This is due to the fact that the dynamically cold 
disc component dominates the central VPs, and to the fact 
that, at larger radii, the second moment of the VP becomes 
large due to the increasing separation of the disc- and bulge- 
components of the VPs. It is remarkable that in most galax- 
ies with a counter-rotating core, the central velocity dis- 
persion strongly increases towards the centre, e.g., IC 1459 
(van der Marel & Franx 1993), NGC 3608 (Jedrzejewski & 
Schechter 1989). More extensive modeling of these systems, 
possibly with a nuclear BH, is clearly required. 

For cases with x ?J 2, the only kinematic signature of 
the counter-rotating disc are strongly skewed VPs. However, 
the interpretation of these as being due to a counter-rotating 
disc is ambiguous. 

6.5 Deriving the central mass density 

The mean streaming of the disc stars follows from the odd 
part of the DF by means of equation (4.5), and substituting 
pd and fg for p and f respectively. This mean streaming 
can be compared to the circular velocity, which is given by 



p(ra 2 )m 2 dm 
y/R 2 ~m 2 (l~ q 2 )' 



(6.6) 



As shown in Section 6.1, the disc rotates with nearly the 
circular velocity (v^, of the disc stars > 0.95u c irc)- Fitting a 
double Gaussian to the observed VPs reveals a hot, mildly 
rotating component, and a rapidly rotating, cold compo- 
nent. The mean rotation of the latter component is an ex- 
cellent measure of the mean streaming of the nuclear disc, 
and therefore of the circular velocity. In turn, this can be 
used to constrain the central mass density. 

We investigated the influence of seeing convolution on 
the rotation velocity of the disc as measured by fitting a dou- 
ble Gaussian to the VPs. The results are shown in Figure 14 
for model 3.1 (a = —1.0, A = 1.0). We plot cod = fdisc/^circ 
as function of the scaling parameter \ f° r three different 
radii: R = 0.5, 1.0, and 2.0 times the scale length of the 



nuclear disc. Seeing convolution decreases uJd considerably, 
especially at R < 2R d . For PSFs whose FWHM exceeds 
2Rd (i.e., \ > 2), the VPs no longer reveal two separate 
components, and the double Gaussian fit becomes meaning- 
less. The fact that u>d ~ 0.8 at R = 0.5Rd for x = (i.e., 
for the unconvolved VPs), is due to the fact that at small 
radii the cusped bulge dominates the surface brightness, and 
therefore these central VPs do not clearly reveal two compo- 
nents, making the interpretation of the double Gaussian fit 
in terms of a disc- and bulge-component somewhat ambigu- 
ous. This is not the case for a non-cusped bulge (a = 0.0), 
where, as shown in Figure 9, the disc dominates the central 
VPs, and has uj d w 0.95. 

7 THE INFLUENCE OF A NUCLEAR BLACK 
HOLE 

Many, if not all, galaxies might harbour massive BHs in 
their nuclei. In order to study their kinematic signature in 
the presence of a nuclear disc, we have constructed several 
models with central BH (see Table 1). A massive BH in the 
nucleus of a galaxy strongly influences the central dynam- 
ics. Early studies based on spherical models showed that 
adiabatic growth of the BH results in a power-law cusp in 
the stellar density profile with a logarithmic slope of central 
surface brightness in the range from —1/2 to —5/4. Further- 
more, hydrostatic equilibrium requires the RMS velocity of 
the stars surrounding the BH to have a r~ 1//2 -cusp (Bahcall 
& Wolff 1976; Young 1980; Quinlan et al. 1995). 

QZMH showed that oblate (a, /3)-models with a BH are 
physical (i.e., have non-negative DF) when a < —0.5. Under 
this condition a nuclear disc can be added, while maintaining 
a positive DF. At sufficiently small radii the potential is 
dominated by the BH, and an explicit expression for p(\&, R) 
can be obtained. Therefore, an asymptotic expression for 
f e (E, L z ) can be calculated in terms of elementary functions 
(see Appendix A) . 

We have added BHs of both M B h = 0.1M d i sc and 
Mbh = -Mdisc to the model with a — —1.0, both for the 
model with and without nuclear disc (see Table 1 for the 
parameters) . The real moments of the VPs after seeing con- 
volution are shown in Figure 15. With this parameterization 
of the VPs the signature of the BH with mass 0.00303Mb u i gc 
is limited to a small increase of the velocity dispersion pro- 
vided that x ^ 0-5- The signature of a BH that is ten 
times more massive, is much more pronounced. Not only is 
a strong central increase of velocity dispersion visible (even 
for x = l)j but in addition, the mean rotation is consider- 
ably larger than in the case without BH. There is no major 
difference of BH-signature between the cases with and with- 
out nuclear disc, as judged from the four velocity profile 
moments presented here. This is due to the fact that the 
nuclear disc only adds a small fraction of the light to the 
VPs, especially at small radii, where the cusp dominates. 

We have also fitted double Gaussians to the seeing con- 
volved VPs. Figure 16 shows the parameter u>d as function of 
radius in arcseconds for three different values of x f° r models 
3.1 (solid lines), 3.3 (dashed lines), and 3.5 (dotted lines). 
The main result is that for the model with Mbh = Mdiac 
(model 3.5) the rotation velocity of the nuclear disc as mea- 
sured from the double Gaussian fit remains over 90% of the 
circular velocity down to 0.5Rd as long as x iS 1- The fact 
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that one can actually measure this is due to the fact that the 
circular velocity in the centre is so large that the small part 
of the VP that is due to the nuclear disc, clearly shows up 
in the wing of the bulge part. In the case of the less massive 
BH (model 3.3), the effect of the BH on the circular velocity 
is so small that the disc part of the VP remains hidden in 
the dominating bulge part, so that the double Gauss fit is 
rather ambiguous, and can no longer be interpreted as due 
to a disc- and bulge-component. 

Since the circular velocity at 0.25i?d for the case with 
Mbh = Afdisc is 4.2 times larger than for the case with- 
out BH, measuring i> c i rc so close to the centre provides an 
excellent measure of the central mass-density and will im- 
mediately reveal the presence of the BH, without the am- 
biguity with respect to velocity anisotropies that hampers 
the interpretation of an observed, central increase of veloc- 
ity dispersion. Therefore, nuclear discs can be used to put 
unambiguous constraints on the presence of nuclear BHs in 
these galaxies. 

8 SUMMARY AND DISCUSSION 

We have constructed axisymmetric, two-integral models of 
elliptical galaxies harbouring small, nuclear discs. We used 
the HQ method to calculate the even part of the phase-space 
distribution function (DF) for models that consist of an (a, 
/3)-spheroid (the bulge) which contains a highly flattened 
exponential spheroid (the disc). In addition, we considered 
the effects of a massive nuclear BH. The bulge has a central 
power- law density cusp (p oc r a ). The disc is thick, and 
has an exponential surface brightness profile. Although we 
concentrated on nuclear discs, the models presented here can 
equally well be used to describe, e.g., SO galaxies, by simply 
increasing the horizontal scale length of the disc. 

For models with only a moderate cusp (i.e., a > —0.5), 
there is a region in the plane of disc-to-bulge ratio A versus 
disc scale-length Rd where the even part of the DF is not 
strictly positive. For the cusped models investigated here 
(a = —0.5, —1.0, —1.5) no such region in (A, Rd) parameter 
space exists, so all these models are physical. We have con- 
sidered different odd parts of the DF, which lead to models 
with a large variety of streaming motions. Upon defining 
the odd part of disc and bulge separately, we obtain disc- 
and bulge-components that are dynamically decoupled. As 
an example of such systems, we have constructed models 
with counter-rotating nuclear discs. 

Our models are locally stable against axisymmetric per- 
turbations according to the Toomre criterion. The nuclear 
discs are stable as long as they are not too flattened or com- 
pact. The presence of a central cusp strongly increases the 
stability of the nuclear disc. Our two-integral nuclear discs 
are also stable against bar formation. Steeper cusps of the 
bulge ensure better stability. The central surface brightness 
of the disc with disc-to-bulge ratio such that the disc is only 
just locally stable over its entire extension depends on the 
value of the disc scale-length R d . The relation is remarkably 
similar to the observed fio-Rd relation for a large sample of 
spirals, SOs, and discy ellipticals. If stability arguments are 
indeed responsible for this relation, this suggests that discs 
build up their mass until they become marginally stable. 

We have investigated the effects of seeing convolution 
on the kinematic observables of nuclear discs. We consid- 



ered PSFs with FWHM in the range of 0.5R d up to 5R d . 
If i?d = 0.2", as observed in a number of Virgo ellipticals, 
this corresponds to the range from 0.1" (i.e., HST resolu- 
tion) up to 1" (i.e., typical ground-based resolution). For 
X = FWHM/i?d > 3 the nuclear disc becomes kinematically 
undetectable. This is similar to the value of \ f° r which 
the disc becomes photometrically undetectable. One does, 
however, observe a slightly larger central velocity dispersion 
than would be the case without nuclear disc. This is due 
to seeing convolution of the rotation curve of the nuclear 
disc. It has been argued by a number of authors that this 
effect might mimic the presence of a central BH. However, 
we find that the effect is rather small, not exceeding 10%, 
and conclude that nuclear discs can not mimic the presence 
of a massive BH. 

In the case of counter-rotating nuclear discs the counter- 
rotation is only unambiguously detectable from the VPs if 
the seeing FWHM is smaller than ~ 2R d (for A = 1). For 
such small seeing FWHM, the central line-of-sight velocity 
profile is dominated by the disc, and therefore one observes 
in addition to the counter-rotation a central decrease in ve- 
locity dispersion. It is surprising that in most cases where 
counter-rotation is observed, one finds an additional strong, 
central increase in velocity dispersion. Although we note 
that this might be an indication for a nuclear BH, in which 
case the high central dispersion is due to seeing convolution 
of the Keplerian rotation curve, further dynamical modeling 
is required to confirm this. 

When fitting a double Gaussian to the VPs one ex- 
tracts information of the kinematics of both components. 
Although the assumption that both components have Gaus- 
sian VPs is clearly simplified, it gives a first order estimate 
of the rotation velocities of the disc. The disc rotates with 
almost the circular velocity (i>di sc ~ 0.95 f c irc)- The rela- 
tive contribution of disc- and bulge light to the VP depends 
strongly on radius, and on the parameters of the model. 
All constructed models have disc-to-bulge ratio chosen so 
that at R = Rd the disc contributes 50% to the light of the 
VP (i.e., A = 1). Therefore, the relative disc contribution 
to the VPs at R < Rd decreases strongly with increasing 
cusp steepness of the bulge. For shallow cusps, the disc still 
contributes a significant amount of light to the central VPs 
and a central decrease in velocity dispersion is observed. 
For more steeply cusped bulges, the disc becomes only de- 
tectable at somewhat larger radii. 

The fact that the nuclear disc rotates with almost the 
circular velocity has important implications. Measuring the 
mean rotation of the nuclear disc allows an accurate determi- 
nation of the central density (or central mass-to-light ratio) 
of its host galaxy. We have tried to "measure" the rotation 
of the nuclear disc by fitting a double Gaussian to the seeing 
convolved VPs. We compared these results to the mean ro- 
tation t>disc as calculated from equation (4.5). For R > 2Rd 
an( l X ^ 1 we found the mean of the rapidly rotating Gaus- 
sian to be in good agreement (within a few percent) with 
fdisc- There VPs can therefore be used to put strong con- 
straints on the central mass-to-light ratio of the host galaxy. 
At smaller radii the cusp dominates and the seeing convolu- 
tion has important influences. Therefore, these central VPs 
do not clearly show two distinct components and the double 
Gaussian fit is somewhat arbitrary. The fit can no longer be 
interpreted as representing the disc- and bulge component. 
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We have also included a massive BH in the centre of 
these models to investigate whether nuclear discs can be 
used to test for the presence of such BHs in the same way 
as ionized gas discs can. In both the cases with and without 
nuclear disc a strong central increase in velocity dispersion is 
observed (for a BH that is massive enough). The observed 
rotation velocities, as measured from the first moment of 
the VP, increases. We have fitted double Gaussians to the 
VPs, to see if we can measure the rotation velocity of the 
nuclear disc. Since the disc rotates with almost the circular 
velocity, such measurements provide an accurate estimate 
of the central mass density of the galaxy. We find that 
for a BH whose mass equals that of the nuclear disc, the 
circular velocity becomes so large close to the centre, that 
even the VPs at only 0.5 ^ from the centre clearly reveal 
two components. The mean velocity of the rapid rotating 
component is an excellent indicator of the circular velocity. 
Since at these small radii, v c i IC is considerably larger than 
for the case without BH, the observed rotation velocity of 
the nuclear disc provides strong evidence for the presence 
of the BH. Most importantly, interpretation of these rota- 
tion velocities is far less ambiguous than interpretation of a 
central increase of velocity dispersion. 

Although the stellar dynamics of the disc and the bulge 
is likely to be influenced by a third integral of motion, we 
suspect that the results discussed here, based upon two- 
integral models, are not strongly influenced by this oversim- 
plification. It is well known that when two-integral models 
require the presence of a massive BH to fit a central increase 
of velocity dispersion, one always has the possibility that a 
three-integral model without BH might equally well fit the 
data. However, if one observes a nuclear disc whose rotation 
velocity exceeds the circular velocity that would correspond 
to a constant mass-to-light ratio model, invoking a third in- 
tegral of motion can not alter the conclusion that there has 
to be a central increase in mass-to-light ratio. 

The high spatial-resolution spectra obtainable with 
HST can be used to accurately measure the rotation curve of 
the nuclear discs, giving an excellent measure of the central 
density of the galaxy. Although the spectral resolution of 
the Faint Object Spectrograph (FOS) aboard the HST may 
be too poor to resolve the VPs in enough detail to allow 
for Udisc to be measured, observations with the Space Tele- 
scope Imaging Spectrograph (STIS), to be installed during 
the next servicing mission, will allow the VPs to be mea- 
sured with high enough both spatial, and spectral resolu- 
tion to put stringent constraints on the presence of possible 
massive BHs in the nuclei of these galaxies. 



ACKNOWLEDGEMENTS 

It is a pleasure to thank Richard Arnold, Marcella Carollo, 
Eric Emsellem, and Eddie Qian for useful discussions. The 
authors are indebted to Jerry Sellwood for helpful advice 
concerning the stability of the nuclear discs. FCvdB was 
supported by the Netherlands Foundation for Astronomi- 
cal Research (ASTRON), #782-373-055, with financial aid 
from the Netherlands Organization for Scientific Research 
(NWO). 



REFERENCES 

Aoki T.E., Hiromoto N., Takami H., Okamura S., 1991, PASJ, 
43, 755 

Athanassoula A., Sellwood J.A., 1986, MNRAS, 221, 213 
Bahcall J.N., Wolff R.A., 1976, ApJ, 209, 214 
Bender R., 1988, A&A, 202, L5 

Bender R., 1990, Dynamics and Interactions of Galaxies, ed. R. 

Wielen (Berlin: Springer Verlag), p. 232 
Binney J. J., Mamon G.A., 1982, MNRAS, 200, 361 
Burkhead M.S., 1986, AJ, 91, 777 

Capaccioli M., Held E.V., Nieto J.-L., 1987, AJ, 94, 1519 
Capaccioli M., Caon N., Rampazzo R., 1990, Dynamics and 

Interactions of Galaxies, ed. R. Wielen (Berlin: Springer 

Verlag), p. 279 

Carollo CM., Franx M., Illingworth G.D., Forbes D.A., 1996, 

ApJ, submitted 
Carter D., 1987, ApJ, 312, 514 

Chandrasekhar, S., 1969, Ellipsoidal Figures of Equilibrium 

(New Haven, Yale University Press) 
Cuddeford P., 1993, MNRAS, 262, 1076 
Dehnen W., Gerhard O.E., 1994, MNRAS, 268, 1019 
Dehnen W., 1995, MNRAS, 274, 919 
Dejonghe H., 1986, Phys. Rep., 133, 217 
de Zeeuw P.T., 1985, MNRAS, 216, 273 
de Zeeuw P.T., Franx M., 1991, ARA&A, 29, 239 
Emsellem E., Monnct G., Bacon R., Nicto J.-L., 1994, A&A, 

285, 739 

Ferrarese L., van den Bosch F.C., Ford H.C., Jaflfe W., 

O'Conncll R.W., 1994, AJ, 108, 1598 
Forbes D.A., 1994, AJ, 107, 2017 

Forbes D.A., Franx M., Illingworth G.D., 1995, AJ, 109, 1988 

Ford H.C., et al., 1994, ApJ, 435, L27 

Freeman K.C., 1970, ApJ, 160, 811 

Gerhard O.E., 1993, MNRAS, 265, 213 

Harms R.J., et al., 1994, ApJ, 435, L35 

Hunter C, 1977, AJ, 82, 271 

Hunter C, Qian E., 1993, MNRAS, 262, 401 (HQ) 
Jacoby G.H., Ciardullo R., Ford H.C., 1990, ApJ, 356, 332 
Jaflfe W., Ford H.C., Ferrarese L., van den Bosch F.C., 

O'Conncll R.W., 1996, ApJ, 460, 214 
Jedrzejewski R.I., 1987, MNRAS, 226, 747 
Jedrzejewski R.I., Schechter P.L., 1989, Dynamics of Dense 

Stellar Systems, ed. D.R. Merritt, (Cambridge: Cambridge 

University Press), p. 25 
Kent S., 1985, ApJS, 59, 115 

Kent S.M., Dame T.M., Fazio G., 1991, ApJ, 378, 131 
Kormendy J., 1988, ApJ, 335, 40 

Kormendy J., Bender R., Richstone D., Ajhar E.A., Dressier A., 
Faber S.M., Gebhardt K., Grillmair C, Laucr T.R., 
Tremainc S., 1996, ApJ, 459, L57 

Kormendy J., Richstone D., 1995, ARA&A, 33, 581 

Kuijken K., 1995, ApJ, 446, 194 

Kuijken K., Gilmore G., 1989, MNRAS, 239, 571 

Laucr T.R., 1985, MNRAS, 216, 429 

Laucr T.R., et al., 1995, AJ, 110, 2622 

Lynden-Bell D., 1962, MNRAS, 124, 1 

Magorrian J., 1995, MNRAS, 277, 1185 

Michard R., 1984, A&A, 140, L39 

Miyoshi M., Moran J., Hcrrnstein J., Greenhill L., Nakai N., 
Diamond P., Inoue M., 1995, Nature, 373, 127 

Nicto J.-L., Poulain P., Davoust E., Rosenblatt P., 1991, A&AS, 
88, 559 

Ostrikcr J. P., Peebles P.J.E., 1973, ApJ, 186, 467 

Qian E., de Zeeuw P.T., van der Marel R.P., Hunter C, 1995, 

MNRAS, 274, 602 (QZMH) 
Quinlan G.D., Hernquist L., Sigurdsson S., 1995, ApJ, 440, 554 
Rix H.-W., White S.D.M., 1990, ApJ, 362, 52 



Models of Elliptical Galaxies with Nuclear Discs 13 



Rix H.-W., White S.D.M., 1992, MNRAS, 254, 389 
Sargent W.L.W., Young P.J., Bokscnberg A., Shortridgc K., 

Lynds C.R., Hartwick F.D.A., 1978, ApJ, 221, 731 
Scorza C, Bender R., 1990, A&A, 235, 49 
Scorza C, Bender R., 1995, A&A, 293, 20 
Seifert W., 1990, PhD Thesis, Univcrsitat Heidelberg 
Shu F.H., 1968, PhD Thesis, Harvard University 
Toomrc A., 1964, ApJ, 139, 1217 

Toomrc A., 1981, The Structure and Evolution of Normal 

Galaxies, eds S.M. Fall, D. Lyndcn-Bcll (Cambridge: 

Cambridge University Press), p. Ill 
van den Bergh S., 1989, PASP, 101, 1072 
van den Bosch F.C., Ferrarese L., Jaffe W., Ford H.C., 

O'Conncll R.W., 1994, AJ, 108, 1579 
van den Bosch F.C., van der Marel R.P., 1995, MNRAS, 274, 

884 

van der Marel R.P., 1991, MNRAS, 253, 710 

van der Marel R.P., 1994, MNRAS, 270, 271 

van der Marel R.P., Franx M., 1993, ApJ, 407, 525 

van der Marel R.P., Evans N.W., Rix H.-W., White S.D.M., de 

Zeeuw P.T., 1994, MNRAS, 271, 99 
Vandervoort P.O., 1970, ApJ, 161, 87 
Young P.J., 1980, ApJ, 242, 1232 

Young P.J., Westphal J. A., Kristian J., Wilson CP., Landauer 
F.P., 1978, ApJ, 221, 721 



where 

ff(£d) =- — t- + 



(i-£ d ) 2 



arctan 



(A5) 



Here £ d = (1 - qj)r/ 2 , and $ d = GMbh / Rd ■ For the radial, 
■q = orbits g(£ d ) = 1- 



APPENDIX: ASYMPTOTIC APPROXIMATION 
OF EVEN PART OF DF CLOSE TO BH 

At sufficiently small radii the density of the composite 
disc/bulge system can be approximated by 

"■>(f) a -"Wi£)' (A1) 

where mb and nid are the constant density spheroids of bulge 
and disc respectively. As long as a > —2 and a + 2/3 < —2, 
the central potential can be approximated by 

•<*"- TO + ,! - (A2 » 

Here is the finite, central stellar potential given by the 
sum of ipo,b (equation 2.9) and ipo,d (equation 2.18). Solv- 
ing z 2 from equation (6.7) and substitution in equation 
(Al) yields an explicit expression of p(\]/,ii). Using the 
alternative form of the contour integral given by equation 
(3.3) in HQ, and integrating along a circular contour pa- 
rameterized by angular parameter 6 (i.e., * = E + Ee 10 , 
~tv < 9 < 7r), an analytic solution for the approximating 
ME, Lz) = f e (E, Lz) + f*{E, L z ) is found. 
For a = — 1 and a + 2/3 = —4 we find 

/ \ - 1 /' 2 



where £ b = (1 — ql)r) 2 , and $ b = GMsii/Rb (see also Dehnen 
and Gerhard 1994, and QZMH). Similarly, for the disc part 
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